### Setup R ###
###############

### Load packages
require(foreign)
require(sandwich)
require(lmtest)
require(MASS)
require(Hmisc)
require(stargazer)
require(sciplot)
require(lme4)

### Clear space
rm(list = ls())

### Load function to incorporate uncertainty (with robust standard errors)
milm <- function(fml, midata){
  xx <- terms(as.formula(fml))
  lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
  vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcovHC(tmp, type = "HC")))
    vcovs[[i]] <- vcovHC(tmp, type = "HC")
  }
  par.est <- apply(lms, 1, mean)
  se.within <- apply(ses, 1, mean)
  se.between <- apply(lms, 1, var)
  se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
  list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)   
}

### Setup Data ###
##################

### Load data. Data primarily from Fariss (2014), Keith et al (2009), and Keith (2012). See Appendix A for data source details and descriptions.
dat <- read.dta("primary.dta")

### Recode NAs
dat$DEMOC[dat$DEMOC==-88] <- NA
dat$DEMOC[dat$DEMOC==-77] <- NA
dat$DEMOC[dat$DEMOC==-66] <- NA
dat$FHStatus[dat$FHStatus==-9] <- NA

### Create democracy dummy variables
dat$POLDEMO <- as.numeric(0)
dat$POLDEMO[dat$DEMOC>5] <- as.numeric(1)
dat$POLDEMO[dat$DEMOC==-88] <- NA
dat$POLDEMO[dat$DEMOC==-77] <- NA
dat$POLDEMO[dat$DEMOC==-66] <- NA
summary(dat$POLDEMO)
dat$FHDEMO <- as.numeric(0)
dat$FHDEMO[dat$FHStatus>1] <- as.numeric(1)
dat$FHDEMO[dat$FHStatus==-9] <- NA
summary(dat$FHDEMO)
dat$GWFDEM <- as.numeric(0)
dat$GWFDEM[dat$gwf_nonautocracy=="democracy"] <- as.numeric(1)
summary(dat$GWFDEM)

### Create logged population variable
dat$logpop <- NA
dat$logpop <- log(dat$POPTOTL)

attach(dat)
datfh <- na.omit(data.frame(latentmean,lag_latentmean,LJI,COWCWAR2,InternationalWar,democracy,MIL2,LEFT,logpop,POPGROWTH, GDPPCPPPCURRENT,GDPPCGROWTH,postsd,lag_latentsd,latentsd,COW,BRIT,uds_mean,uds_sd,POLDEMO,FHDEMO,GWFDEM,YEAR,DEMOC,KeithJUDIND,LJIwoCIM,LJIwoCIMpostsd,YEAR,COW))
 dim(datfh)

### Analyze Data ###
####################

### Estimate model with point estimates from latent variables
lmdd.out <- lm(latentmean ~ lag_latentmean + LJI + 
COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH, data=datfh)
summary(lmdd.out)
lmdd.out$robustse <- vcovHC(lmdd.out)
a <- coeftest(lmdd.out,lmdd.out$robustse)
a

### Plot estimates
var.names <- c("De Facto Judicial Independence \n (Latent Measure)","Civil War","International War","Democracy","Military Control","State-Socialist Regime",
               "British Col. Exper.","Economic Development","Economic Growth","Logged Population","Population Growth")
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.1, .1),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA)
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- coef(lmdd.out)[c(-1,-2)]
se <- sqrt(diag(vcovHC(lmdd.out)))[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#34495e")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#2980b9")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")
    
### Setup data to run models where we account for uncertainty
set.seed(12261991)
newdata <- list()
for(i in 1:1000){
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$latentmean, sd=datfh$latentsd)
  newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_latentmean, sd=datfh$lag_latentsd)
  newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJI, sd=datfh$postsd)
    newdata[[i]]$ljiwocim <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJIwoCIM, sd=datfh$LJIwoCIMpostsd)
  newdata[[i]]$uds <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$uds_mean, sd=datfh$uds_sd)
}

### Check with uncertainty in DV and lagged DV
FML <- "draw ~ drawlag + LJI + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3002
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Plot estimates
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.75, .75),
     ylim = c(1, length(var.names)),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- model$beta[c(-1,-2)]
se <- model$SE[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#2c3e50")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3) 
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#d35400")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")

### Check with uncertainty in DV, lagged DV, and IV
FML2 <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model2 <- milm(fml=FML2, midata=newdata)
df <- 3002
tstats2<-data.frame(model2$terms,model2$beta,model2$SE,model2$beta/model2$SE,2*pt(-abs(model2$beta/model2$SE),df=df))
colnames(tstats2) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats2

### Plot estimates
plot(NULL,
     xlim = c(-.75, .75),
     ylim = c(.8, length(var.names)),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- model2$beta[c(-1,-2)]
se <- model2$SE[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#2c3e50")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#27ae60")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")

### Plot estimates from all 3 models
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.5, .5),
     ylim = c(0.9, 3),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 3) 
# without uncertainty
lines(c(0.05691, -0.02717), c(3, 3), lwd = 4, col="#95a5a6")
lines(c(0.05005, -0.0203), c(3, 3), lwd = 10, col="#2980b9")
points(0.01487271, 3, pch = 19, cex = 1, col="white")
points(0.01487271, 3, pch = 19, cex = .5, col="black")
text(0.01487271, 3, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.01487271, 3, c("(Model 1)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in DV, and lagged DV
lines(c(0.3532, 0.09939), c(2, 2), lwd = 4, col="#95a5a6")
lines(c(0.3325, 0.1201), c(2, 2), lwd = 10, col="#d35400")
points(0.226284542, 2, pch = 19, cex = 1, col="white")
points(0.226284542, 2, pch = 19, cex = .5, col="black")
text(0.226284542, 2, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.226284542, 2, c("(Model 2 - Uncertainty in DV, and Lagged DV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in DV, lagged DV, and IV
lines(c(0.3224, 0.08377), c(1, 1), lwd = 4, col="#95a5a6")
lines(c(0.3029, 0.1032), c(1, 1), lwd = 10, col="#27ae60")
points(0.20309479, 1, pch = 19, cex = 1, col="white")
points(0.20309479, 1, pch = 19, cex = .5, col="black")
text(0.20309479, 1, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.20309479, 1, c("(Model 3 - Uncertainty in DV, Lagged DV, and IV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 

### Cross Validation ###
########################

### Declare formulas and estimate OLS models
FORMULA <- "draw ~ drawlag"

model <- milm(fml=FORMULA, midata=newdata)
df <- 3113
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA2 <- "draw ~ drawlag + ljidraw"

model <- milm(fml=FORMULA2, midata=newdata)
df <- 3112
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA3 <- "draw ~ drawlag + COWCWAR2 + InternationalWar"

model <- milm(fml=FORMULA3, midata=newdata)
df <- 3111
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA4 <- "draw ~ drawlag + democracy + MIL2 + LEFT + BRIT"

model <- milm(fml=FORMULA4, midata=newdata)
df <- 3009
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA5 <- "draw ~ drawlag + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"

model <- milm(fml=FORMULA5, midata=newdata)
df <- 3009
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA6 <- "draw ~ drawlag + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"

model <- milm(fml=FORMULA6, midata=newdata)
df <- 3003
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

FORMULA7 <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"

model <- milm(fml=FORMULA7, midata=newdata)
df <- 3002
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Declare functions
loss.rmse <- function(fit, df) cbind(as.numeric(predict(fit, newdata = df)), df$draw)

model.add <- function(df) lm(FORMULA, data=df)
model.add2 <- function(df) lm(FORMULA2, data=df)
model.add3 <- function(df) lm(FORMULA3, data=df)
model.add4 <- function(df) lm(FORMULA4, data=df)
model.add5 <- function(df) lm(FORMULA5, data=df)
model.add6 <- function(df) lm(FORMULA6, data=df)
model.add7 <- function(df) lm(FORMULA7, data=df)

validate.cv <- function(df, folds, resamples, model, loss) {
  
  lapply(1:resamples, function(x) {
    
    df$folds <- sample(rep(1:folds, length.out = nrow(df)))
    
    lapply(1:folds, function(test) {
      
      fit <- model(df[df$folds != test, ])
      
      loss(fit, df[df$folds == test, ])
      
    })})}

### Set values
SIMS <- 1000
cv1 <- NA
cv.add <- NA
cv2 <- NA
cv.add2 <- NA
cv3 <- NA
cv.add3 <- NA
cv4 <- NA
cv.add4 <- NA
cv5 <- NA
cv.add5 <- NA
cv6 <- NA
cv.add6 <- NA
cv7 <- NA
cv.add7 <- NA
newdata <- list()

for(i in 1:SIMS){
  
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$latentmean, sd=datfh$latentsd)
  newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_latentmean, sd=datfh$lag_latentsd)
  newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJI, sd=datfh$postsd)
      newdata[[i]]$ljiwocim <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJIwoCIM, sd=datfh$LJIwoCIMpostsd)
  
  cv.add <- validate.cv(newdata[[i]], 10, 1, model.add, loss.rmse)
  cv.add2 <- validate.cv(newdata[[i]], 10, 1, model.add2, loss.rmse)
  cv.add3 <- validate.cv(newdata[[i]], 10, 1, model.add3, loss.rmse)
  cv.add4 <- validate.cv(newdata[[i]], 10, 1, model.add4, loss.rmse)  
  cv.add5 <- validate.cv(newdata[[i]], 10, 1, model.add5, loss.rmse)  
  cv.add6 <- validate.cv(newdata[[i]], 10, 1, model.add6, loss.rmse)  
  cv.add7 <- validate.cv(newdata[[i]], 10, 1, model.add7, loss.rmse)  
  
  value <- rbind(cv.add[[1]][[1]],cv.add[[1]][[2]],cv.add[[1]][[3]],cv.add[[1]][[4]],cv.add[[1]][[5]],cv.add[[1]][[6]],cv.add[[1]][[7]],cv.add[[1]][[8]],cv.add[[1]][[9]],cv.add[[1]][[10]])
  cv1[i] <- sqrt(mean((value[,1] - value[,2])^2))
  value2 <- rbind(cv.add2[[1]][[1]],cv.add2[[1]][[2]],cv.add2[[1]][[3]],cv.add2[[1]][[4]],cv.add2[[1]][[5]],cv.add2[[1]][[6]],cv.add2[[1]][[7]],cv.add2[[1]][[8]],cv.add2[[1]][[9]],cv.add2[[1]][[10]])
  cv2[i] <- sqrt(mean((value2[,1] - value2[,2])^2))
  value3 <- rbind(cv.add3[[1]][[1]],cv.add3[[1]][[2]],cv.add3[[1]][[3]],cv.add3[[1]][[4]],cv.add3[[1]][[5]],cv.add3[[1]][[6]],cv.add3[[1]][[7]],cv.add3[[1]][[8]],cv.add3[[1]][[9]],cv.add3[[1]][[10]])
  cv3[i] <- sqrt(mean((value3[,1] - value3[,2])^2))
  value4 <- rbind(cv.add4[[1]][[1]],cv.add4[[1]][[2]],cv.add4[[1]][[3]],cv.add4[[1]][[4]],cv.add4[[1]][[5]],cv.add4[[1]][[6]],cv.add4[[1]][[7]],cv.add4[[1]][[8]],cv.add4[[1]][[9]],cv.add4[[1]][[10]])
  cv4[i] <- sqrt(mean((value4[,1] - value4[,2])^2))
  value5 <- rbind(cv.add5[[1]][[1]],cv.add5[[1]][[2]],cv.add5[[1]][[3]],cv.add5[[1]][[4]],cv.add5[[1]][[5]],cv.add5[[1]][[6]],cv.add5[[1]][[7]],cv.add5[[1]][[8]],cv.add5[[1]][[9]],cv.add5[[1]][[10]])
  cv5[i] <- sqrt(mean((value5[,1] - value5[,2])^2))
  value6 <- rbind(cv.add6[[1]][[1]],cv.add6[[1]][[2]],cv.add6[[1]][[3]],cv.add6[[1]][[4]],cv.add6[[1]][[5]],cv.add6[[1]][[6]],cv.add6[[1]][[7]],cv.add6[[1]][[8]],cv.add6[[1]][[9]],cv.add6[[1]][[10]])
  cv6[i] <- sqrt(mean((value6[,1] - value6[,2])^2))
  value7 <- rbind(cv.add7[[1]][[1]],cv.add7[[1]][[2]],cv.add7[[1]][[3]],cv.add7[[1]][[4]],cv.add7[[1]][[5]],cv.add7[[1]][[6]],cv.add7[[1]][[7]],cv.add7[[1]][[8]],cv.add7[[1]][[9]],cv.add7[[1]][[10]])
  cv7[i] <- sqrt(mean((value7[,1] - value7[,2])^2))
}

y <- c(((cv2-cv1)/cv1),((cv3-cv1)/cv1),((cv4-cv1)/cv1),((cv5-cv1)/cv1),((cv6-cv1)/cv1),((cv7-cv1)/cv1))
a <- c(rep(c("(2) De Facto Judicial \n Independence \n Only"),1000),rep(c("(3) Threats \n Only \n"),1000),rep(c("(4) Regime Type \n Only \n "),1000),
       rep(c("(5) Socioeconomic \n Only \n"),1000),rep(c("(6) All Controls \n \n"),1000),rep(c("(6) Full Model \n \n"),1000))
d <- cbind(as.numeric(y),a)

### Plot reduction in mean square error across models
par(
  oma = c(0,0,0,0),
  mar = c(5,4,2,1)
)
bargraph.CI(d[,2],as.numeric(d[,1]), main="", col=c("#d35400","#cccccc","#cccccc","#cccccc","#cccccc","#d35400"), font.main=1,
            border=NA, family="Helvetica", ylim=c(min(as.numeric(d[,1])),0),err.width=.1,xlab="Models",ylab="Proportional Reduction in Mean Square Error")

#########################
### Robustness Checks ###
#########################

### Appendix B ###
##################

### Examine correlations between latent judicial independence measure and democracy measures                        
d <- data.frame(x1=datfh$LJI,x2=datfh$democracy,x3=datfh$DEMOC,x4=datfh$POLDEMO,x5=datfh$FHDEMO,x6=datfh$GWFDEM,x7=datfh$uds_mean)
cor(d)

### Appendix C ###
##################

### Check with uncertainty in DV and lagged DV but with modified LJI measure
FML2 <- "draw ~ drawlag + ljiwocim + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model2 <- milm(fml=FML2, midata=newdata)
df <- 3002
tstats2<-data.frame(model2$terms,model2$beta,model2$SE,model2$beta/model2$SE,2*pt(-abs(model2$beta/model2$SE),df=df))
colnames(tstats2) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats2

### Plot estimates
plot(NULL,
     xlim = c(-.75, .75),
     ylim = c(.8, length(var.names)),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
est <- model2$beta[c(-1,-2)]
se <- model2$SE[c(-1,-2)]
for (i in 1:length(est)) {
  lines(c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), c(i, i), lwd = 2, col="#95a5a6")
  lines(c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), c(i, i), lwd = 5, col="#2c3e50")
  points(est[i], i, pch = 19, cex = 1, col="white")
  points(est[i], i, pch = 19, cex = .5, col="black")
  text(est[i], i, var.names[i], xpd = T, cex = .8, pos = 3, offset=0.3)
}
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 2.5) 
lines(c(est[1] + 1.96*se[1], est[1] - 1.96*se[1]), c(1, 1), lwd = 2, col="#95a5a6")
lines(c(est[1] + 1.64*se[1], est[1] - 1.64*se[1]), c(1, 1), lwd = 5, col="#27ae60")
points(est[1], 1, pch = 19, cex = 1, col="white")
points(est[1], 1, pch = 19, cex = .5, col="black")
text(est[1], 0.4, c("(Without CIM)"), xpd = T, cex = .8, pos = 3, offset=0.5) 

### Appendix D ###
##################

### Check with uncertainty in IV
FML <- "latentmean ~ lag_latentmean + ljidraw + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3002
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Plot point estimates from all 4 models
par(
  oma = c(0,0,0,0),
  mar = c(5,2,4,2)
)
plot(NULL,
     xlim = c(-.5, .5),
     ylim = c(0.9, 4),
     axes = F, xlab = NA, ylab = NA) 
abline(v = 0, lty = 3, lwd=2, col = "#cccccc")
axis(side = 1)
mtext(side = 1, "Estimated Coefficients", line = 3) 
# without uncertainty
lines(c(0.05691, -0.02717), c(4, 4), lwd = 4, col="#95a5a6")
lines(c(0.05005, -0.0203), c(4, 4), lwd = 10, col="#2980b9")
points(0.01487271, 4, pch = 19, cex = 1, col="white")
points(0.01487271, 4, pch = 19, cex = .5, col="black")
text(0.01487271, 4, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.01487271, 4, c("(Model 1)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in IV
lines(c(0.05443, -0.02797), c(3, 3), lwd = 4, col="#95a5a6")
lines(c(0.0477, -0.02124), c(3, 3), lwd = 10, col="#CCCCFF")
points(0.013229894, 3, pch = 19, cex = 1, col="white")
points(0.013229894, 3, pch = 19, cex = .5, col="black")
text(0.013229894, 3, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.013229894, 3, c("(Model 2 - Uncertainty in IV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in DV, and lagged DV
lines(c(0.3532, 0.09939), c(2, 2), lwd = 4, col="#95a5a6")
lines(c(0.3325, 0.1201), c(2, 2), lwd = 10, col="#d35400")
points(0.226284542, 2, pch = 19, cex = 1, col="white")
points(0.226284542, 2, pch = 19, cex = .5, col="black")
text(0.226284542, 2, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.226284542, 2, c("(Model 3 - Uncertainty in DV, and Lagged DV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 
# uncertainty in DV, lagged DV, and IV
lines(c(0.3224, 0.08377), c(1, 1), lwd = 4, col="#95a5a6")
lines(c(0.3029, 0.1032), c(1, 1), lwd = 10, col="#27ae60")
points(0.20309479, 1, pch = 19, cex = 1, col="white")
points(0.20309479, 1, pch = 19, cex = .5, col="black")
text(0.20309479, 1, c("De Facto Judicial Independence \n (Latent Measure)"), xpd = T, cex = .8, pos = 3, offset=0.5) 
text(0.20309479, 1, c("(Model 4 - Uncertainty in DV, Lagged DV, and IV)"), xpd = T, cex = .8, pos = 1, offset=0.5) 

### Appendix F ###
##################

### Check with LJI measure that incorporates uncertainty in the DV, lagged DV, and IV using the UDS measure
FML3 <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + uds + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model3 <- milm(fml=FML3, midata=newdata)
df <- 2238
tstats3<-data.frame(model3$terms,model3$beta,model3$SE,model3$beta/model3$SE,2*pt(-abs(model3$beta/model3$SE),df=df))
colnames(tstats3) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats3

#########################
### Unreported Models ###
#########################

### Create new data frame
datfh <- na.omit(data.frame(latentmean,lag_latentmean,LJI,COWCWAR2,InternationalWar,democracy,MIL2,LEFT,logpop,POPGROWTH, GDPPCPPPCURRENT,GDPPCGROWTH,postsd,lag_latentsd,latentsd,COW,BRIT,uds_mean,uds_sd,POLDEMO,FHDEMO,GWFDEM,YEAR,DEMOC,KeithJUDIND,LJIwoCIM,LJIwoCIMpostsd,YEAR,lag_LJI,lag_postsd))
dim(datfh)

### Examine Effect of Lagged Judicial Independence ###
######################################################

### Setup data to run models when we account for uncertainty
set.seed(12261991)
newdata <- list()
for(i in 1:1000){
  newdata[[i]] <- datfh
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$latentmean, sd=datfh$latentsd)
  newdata[[i]]$drawlag <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_latentmean, sd=datfh$lag_latentsd)
  newdata[[i]]$ljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$LJI, sd=datfh$postsd)
    newdata[[i]]$lagljidraw <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$lag_LJI, sd=datfh$lag_postsd)
  newdata[[i]]$uds <- rnorm(cbind(rep(1,nrow(datfh))), mean=datfh$uds_mean, sd=datfh$uds_sd)
}

### Check with uncertainty in all using lagged IV
FML <- "draw ~ drawlag + lagljidraw + COWCWAR2 + InternationalWar + democracy + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3302
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Check with uncertainty in all using POLDEMO
FML <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + POLDEMO + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3302
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Check with uncertainty in all using FHDEMO
FML <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + FHDEMO + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3302
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Check with uncertainty in all using DEMOC
FML <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + DEMOC + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3302
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats

### Check with uncertainty in all using GWFDEM
FML <- "draw ~ drawlag + ljidraw + COWCWAR2 + InternationalWar + GWFDEM + MIL2 + LEFT + BRIT + GDPPCPPPCURRENT + GDPPCGROWTH + logpop + POPGROWTH"
model <- milm(fml=FML, midata=newdata)
df <- 3302
tstats<-data.frame(model$terms,model$beta,model$SE,model$beta/model$SE,2*pt(-abs(model$beta/model$SE),df=df))
colnames(tstats) <- c("Variable","Est. Coef.","SE","t-stat","p")
tstats